# =============================================================================
# MAIN FIGURES 1-6: VISUALIZATIONS FOR MAIN TEXT
# Creates all primary figures for the main text of the paper
# =============================================================================

# --- PACKAGES ----------------------------------------------------------------
library(readr)
library(haven)
library(dplyr)
library(ggplot2)
library(tidyr)

# --- DATA LOADING ------------------------------------------------------------
vdem      <- read_csv("vdem.csv")
miller    <- read_dta("AutocraticRulingPartiesDataset_cy.dta")
estimates <- read_csv("estimates_independent.csv")

# =============================================================================
# FIGURE 1: PREVALENCE OF INSTITUTIONS (1946-2015)
# =============================================================================

# --- DATA PREPARATION --------------------------------------------------------
# V-Dem subset: autocracies only (v2x_regime < 2), years > 1945
vdem_subset <- vdem %>%
  select(COWcode, year, v2lgbicam, v2x_regime, v2juhcind_ord) %>%
  filter(v2x_regime < 2, year > 1945) %>%
  mutate(
    leg = ifelse(v2lgbicam > 0, 1, 0),
    jud = ifelse(v2juhcind_ord > 1, 1, 0)
  )

# Miller dataset preparation
miller_subset <- miller %>%
  rename(COWcode = ccode, year = Year) %>%
  filter(year > 1945)

# Merge datasets
merged_data <- inner_join(vdem_subset, miller_subset, by = c("COWcode", "year"))

# Calculate shares by year
leg_share_by_year   <- merged_data %>% 
  group_by(year) %>% 
  summarise(share_leg = mean(leg, na.rm = TRUE), .groups = "drop")

party_share_by_year <- merged_data %>% 
  group_by(year) %>% 
  summarise(share_party = mean(Party_Regime, na.rm = TRUE), .groups = "drop")

jud_share_by_year   <- merged_data %>% 
  group_by(year) %>% 
  summarise(share_jud = mean(jud, na.rm = TRUE), .groups = "drop")

combined_share <- leg_share_by_year %>%
  left_join(party_share_by_year, by = "year") %>%
  left_join(jud_share_by_year,   by = "year") %>%
  pivot_longer(cols = c(share_leg, share_party, share_jud),
               names_to = "Category", values_to = "Share")

# --- VISUALIZATION -----------------------------------------------------------
p_institutions <- ggplot(combined_share, aes(x = year, y = Share, linetype = Category)) +
  geom_line(color = "black", linewidth = 0.8) +
  scale_x_continuous(breaks = seq(1946, 2015, by = 4)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.1), limits = c(0, 1)) +
  scale_linetype_manual(values = c("share_jud" = "solid", "share_leg" = "twodash", 
                                   "share_party" = "dotted"),
                        labels = c("Judiciary", "Legislature", "Ruling Party")) +
  labs(x = "Year", y = "Share of Autocracies", linetype = "") +
  theme_minimal(base_size = 11) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.line = element_line(color = "black"), legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1))

# --- OUTPUT ------------------------------------------------------------------
ggsave("Figure1.pdf", p_institutions)

# =============================================================================
# FIGURE 2: DENSITY BY REGIME TYPE
# =============================================================================

# --- VISUALIZATION -----------------------------------------------------------
p_density <- ggplot(estimates, aes(x = dyn.estimates, y = after_stat(density), 
                                   fill = factor(v2x_regime))) +
  geom_density(alpha = 0.5) +
  scale_fill_manual(values = c("0" = "grey50", "1" = "grey80"),
                    name = "", labels = c("Closed Autocracies", "Electoral Autocracies")) +
  scale_x_continuous(limits = c(-8, 4), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, 0.8), expand = c(0, 0)) +
  labs(x = expression(Latent~estimate~theta), y = "Density", title = NULL) +
  theme_classic() +
  theme(legend.position = "bottom")

# --- OUTPUT ------------------------------------------------------------------
ggsave("Figure2.pdf", p_density)

# =============================================================================
# FIGURE 3: TREND OVER TIME
# =============================================================================

# --- DATA PREPARATION --------------------------------------------------------
yearly_avg <- estimates %>%
  group_by(year) %>%
  summarise(
    mean_dyn_estimates = mean(dyn.estimates, na.rm = TRUE),
    mean_dyn_up = mean(dyn.up, na.rm = TRUE),
    mean_dyn_lo = mean(dyn.lo, na.rm = TRUE), 
    .groups = "drop"
  )

# --- VISUALIZATION -----------------------------------------------------------
p_mean <- ggplot(yearly_avg, aes(x = year, y = mean_dyn_estimates)) +
  geom_ribbon(aes(ymin = mean_dyn_lo, ymax = mean_dyn_up), alpha = 0.3, color = NA) +
  geom_line(linewidth = 0.75) +
  scale_x_continuous(breaks = seq(1946, 2023, by = 2)) +
  ylim(-1, 0.75) +
  geom_vline(xintercept = 1975, linetype = "dotted", linewidth = 1) +
  labs(x = "Year", y = "Latent Estimate", title = NULL) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# --- OUTPUT ------------------------------------------------------------------
ggsave("Figure3.pdf", p_mean)

# =============================================================================
# FIGURE 4: CROSS-SECTION 2023
# =============================================================================

# --- DATA PREPARATION --------------------------------------------------------
latest <- estimates %>%
  filter(year == 2023, !is.na(country_name), !is.na(dyn.estimates)) %>%
  arrange(dyn.estimates, country_name) %>%
  mutate(country_factor = factor(country_name, levels = unique(country_name), ordered = TRUE))

# --- VISUALIZATION -----------------------------------------------------------
p_2023 <- ggplot(latest, aes(x = dyn.estimates, y = country_factor)) +
  geom_segment(aes(x = dyn.lo, xend = dyn.up, y = country_factor, yend = country_factor),
               linewidth = 0.4, alpha = 0.9, na.rm = TRUE) +
  geom_point(size = 1.5) +
  scale_x_continuous(limits = c(-10, 5), breaks = seq(-10, 5, by = 5)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  labs(x = expression(Latent~estimate~theta), y = "Country", title = NULL) +
  theme_bw() +
  theme(
    panel.border = element_rect(color = "black", linewidth = 1),
    axis.line = element_line(colour = "black"),
    axis.ticks = element_line(colour = "black"),
    axis.text.y = element_text(size = 9)
  )

# --- OUTPUT ------------------------------------------------------------------
ggsave("Figure4.pdf", p_2023)

# =============================================================================
# FIGURE 5: KENYA CASE STUDY
# =============================================================================

# --- DATA PREPARATION --------------------------------------------------------
subset_data_kenya <- estimates %>%
  filter(country_name == "Kenya", year > 1963, year < 2003) %>%
  select(year, country_name, dyn.estimates, dyn.up, dyn.lo)

# --- VISUALIZATION -----------------------------------------------------------
p_kenya <- ggplot(subset_data_kenya, aes(x = year, y = dyn.estimates)) +
  geom_line() +
  geom_ribbon(aes(ymin = dyn.lo, ymax = dyn.up), alpha = 0.2) +
  labs(x = "Year", y = expression(Latent~Estimate~theta), title = "") +
  scale_x_continuous(breaks = seq(1964, 2002, by = 2)) +
  scale_y_continuous(limits = c(-1, 2), breaks = seq(-1, 2, by = 0.5)) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
    plot.title = element_text(hjust = 0.5)
  ) +
  geom_vline(xintercept = 1964, linetype = "dotted", linewidth = 0.5) +
  geom_text(label = "Jomo Kenyatta (1964-1978)", x = 1964.75, y = 1.4, hjust = 0, size = 4) +
  geom_vline(xintercept = 1978, linetype = "dotted", linewidth = 0.5) +
  geom_text(label = "Daniel arap Moi (1978-2002)", x = 1983.5, y = 1.4, hjust = 0, size = 4) +
  geom_vline(xintercept = 2002, linetype = "dotted", linewidth = 0.5)

# --- OUTPUT ------------------------------------------------------------------
ggsave("Figure5.pdf", plot = p_kenya)

# =============================================================================
# FIGURE 6: ZAMBIA CASE STUDY
# =============================================================================

# --- DATA PREPARATION --------------------------------------------------------
subset_data_zambia <- estimates %>%
  filter(country_name == "Zambia", year < 1992) %>%
  select(year, country_name, dyn.estimates, dyn.up, dyn.lo)

# --- VISUALIZATION -----------------------------------------------------------
p_zambia <- ggplot(subset_data_zambia, aes(x = year, y = dyn.estimates)) +
  geom_line() +
  geom_ribbon(aes(ymin = dyn.lo, ymax = dyn.up), alpha = 0.2) +
  labs(x = "Year", y = expression(Latent~Estimate~theta), title = "") +
  scale_x_continuous(breaks = seq(1964, 1991, by = 2)) +
  scale_y_continuous(limits = c(-1, 1), breaks = seq(-1, 1, by = 0.5)) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
    plot.title = element_text(hjust = 0.5)
  ) +
  geom_text(label = "Kenneth Kaunda (1964-1991)", x = 1973.5, y = 0.5, hjust = 0, size = 4)

# --- OUTPUT ------------------------------------------------------------------
ggsave("Figure6.pdf", plot = p_zambia)